import pandas as pd
import numpy as np
import geopandas as gpd
import matplotlib.pyplot as plt
import seaborn as sns
import pysal as py
from matplotlib import colors
from pysal.viz import splot
from esda.moran import Moran_Local
from splot.esda import plot_moran
sns.set_context(context='paper')
%matplotlib inline
from census import Census
from us import states
from shapely.geometry import Point
from pyproj import Proj, transform
from geopy.distance import great_circle
from pysal.explore import esda
from pysal.lib import weights
from numpy.random import seed
import rioxarray  # Surface data manipulation
import xarray  # Surface data manipulation
from pysal.model import spreg
import statsmodels.formula.api as sm
## warnings is a module that allows you to filter warnings
import warnings
## we are going to ignore all warnings (so they won't print)
warnings.filterwarnings("ignore")

C:\Users\86158\anaconda3\envs\gds_py\lib\site-packages\numba\core\decorators.py:262: NumbaDeprecationWarning: numba.generated_jit is deprecated. Please see the documentation at: https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-generated-jit for more information and advice on a suitable replacement.
  warnings.warn(msg, NumbaDeprecationWarning)
C:\Users\86158\anaconda3\envs\gds_py\lib\site-packages\quantecon\lss.py:20: NumbaDeprecationWarning: The 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.
  def simulate_linear_model(A, x0, v, ts_length):
C:\Users\86158\anaconda3\envs\gds_py\lib\site-packages\spaghetti\network.py:40: FutureWarning: The next major release of pysal/spaghetti (2.0.0) will drop support for all ``libpysal.cg`` geometries. This change is a first step in refactoring ``spaghetti`` that is expected to result in dramatically reduced runtimes for network instantiation and operations. Users currently requiring network and point pattern input as ``libpysal.cg`` geometries should prepare for this simply by converting to ``shapely`` geometries.
  warnings.warn(dep_msg, FutureWarning, stacklevel=1)

c = Census("d9c002dc1334c8f6cbea48d3f10a4176cdf89064")

# Retrieving census data for Cook County, Illinois, for the year 2021
_2021_census = c.acs5.state_county_tract(fields=('NAME', 'B25077_001E', 'B25085_001E', 'B25001_001E', 'B25002_003E', 'B25003_002E', 'B25003_003E', 'B06011_001E', 'B05007_001E', 'B02001_001E','B02001_002E'),
                                           state_fips=states.IL.fips,  
                                           county_fips="031",  # Cook County, Illinois
                                           tract="*",
                                           year=2021)
# Retrieving census data for Cook County, Illinois, for the year 2021
_2018_census = c.acs5.state_county_tract(fields=('NAME', 'B25077_001E', 'B25085_001E', 'B25001_001E', 'B25002_003E', 'B25003_002E', 'B25003_003E', 'B06011_001E', 'B05007_001E', 'B02001_001E','B02001_002E'),
                                           state_fips=states.IL.fips,  
                                           county_fips="031",  # Cook County, Illinois
                                           tract="*",
                                           year=2018)

#_2021_census

# Convert the data to a pandas dataframe
acs_df_2021 = pd.DataFrame(_2021_census)
# Create a GEOID column by concatenating the state, county, and tract codes
acs_df_2021["GEOID"] = acs_df_2021["state"] + acs_df_2021["county"] + acs_df_2021["tract"]
# Calculate the percentage of non-White population by subtracting the White population from the total population and dividing by the total population
acs_df_2021["nonwhite_perc"] = (acs_df_2021["B02001_001E"] - acs_df_2021["B02001_002E"]) / acs_df_2021["B02001_001E"]
acs_df_2021["foreign_perc"] = acs_df_2021["B05007_001E"] / acs_df_2021["B02001_001E"]

acs_df_2021["vacant_unit_perc"] = acs_df_2021["B25002_003E"] / acs_df_2021["B25001_001E"]
acs_df_2021["owner_occupied_perc"] = acs_df_2021["B25003_002E"] / acs_df_2021["B25001_001E"]
acs_df_2021["renter_occupied_perc"] = acs_df_2021["B25003_003E"] / acs_df_2021["B25001_001E"]

# Convert the data to a pandas dataframe
acs_df_2018 = pd.DataFrame(_2018_census)
# Create a GEOID column by concatenating the state, county, and tract codes
acs_df_2018["GEOID"] = acs_df_2018["state"] + acs_df_2018["county"] + acs_df_2018["tract"]
# Calculate the percentage of non-White population by subtracting the White population from the total population and dividing by the total population
acs_df_2018["nonwhite_perc"] = (acs_df_2018["B02001_001E"] - acs_df_2018["B02001_002E"]) / acs_df_2018["B02001_001E"]
acs_df_2018["foreign_perc"] = acs_df_2018["B05007_001E"] / acs_df_2018["B02001_001E"]

acs_df_2018["vacant_unit_perc"] = acs_df_2018["B25002_003E"] / acs_df_2018["B25001_001E"]
acs_df_2018["owner_occupied_perc"] = acs_df_2018["B25003_002E"] / acs_df_2018["B25001_001E"]
acs_df_2018["renter_occupied_perc"] = acs_df_2018["B25003_003E"] / acs_df_2018["B25001_001E"]

# Read the shapefile for Cook County
cook_tract = gpd.read_file("https://www2.census.gov/geo/tiger/TIGER2021/TRACT/tl_2021_17_tract.zip")
# Filter the tracts for Cook County only
cook_tract = cook_tract[cook_tract["COUNTYFP"] == "031"]

cook_census_geo_2021 = cook_tract.merge(acs_df_2021, left_on = 'GEOID', right_on = 'GEOID')
cook_census_geo_2021.describe()

cook_census_geo_2018 = cook_tract.merge(acs_df_2018, left_on = 'GEOID', right_on = 'GEOID')
cook_census_geo_2018.describe()

cook_census_geo_2018.head()

file_path = "data/Boundaries - Neighborhoods.geojson"
chicago_boundary = gpd.read_file(file_path)
chicago_boundary = chicago_boundary.to_crs(cook_census_geo_2021.crs)
cook_census_geo_2021_clipped = gpd.clip(cook_census_geo_2021, chicago_boundary)
cook_census_geo_2018_clipped = gpd.clip(cook_census_geo_2018, chicago_boundary)

# Create a figure and axis object
fig, axs = plt.subplots(2, 4, figsize=(20, 8))
# Plot 2021 data
columns_2021 = ['nonwhite_perc', 'vacant_unit_perc', 'foreign_perc', 'renter_occupied_perc']
titles_2021 = ['Non-White Percentage by Census Tract (2021)', 'Vacant Unit Percentage by Census Tract (2021)', 'Foreign-Born Percentage by Census Tract (2021)', 'Renter Occupied Percentage by Census Tract (2021)']
for i in range(len(columns_2021)):
    cook_census_geo_2021_clipped.plot(column=columns_2021[i], cmap='coolwarm', linewidth=0.8, ax=axs[0, i], edgecolor='0.8', legend=True)
    axs[0, i].set_title(titles_2021[i])
# Plot 2018 data
columns_2018 = ['nonwhite_perc', 'vacant_unit_perc', 'foreign_perc', 'renter_occupied_perc']
titles_2018 = ['Non-White Percentage by Census Tract (2018)', 'Vacant Unit Percentage by Census Tract (2018)', 'Foreign-Born Percentage by Census Tract (2018)', 'Renter Occupied Percentage by Census Tract (2018)']
for i in range(len(columns_2018)):
    cook_census_geo_2018_clipped.plot(column=columns_2018[i], cmap='coolwarm', linewidth=0.8, ax=axs[1, i], edgecolor='0.8', legend=True)
    axs[1, i].set_title(titles_2018[i])
plt.tight_layout()
plt.show()

# Create a figure and axis object
fig, ax = plt.subplots(1, 2, figsize=(18, 9))
# Plot the Choropleth map with the filtered data for 2021
cook_census_geo_2021_filtered_clipped.plot(column='B25077_001E', ax=ax[0], legend=True, scheme='Quantiles', cmap='coolwarm')
ax[0].set_title('Median Value of Owner-Occupied Housing Units (2021)', fontsize=16)
# Plot the Choropleth map with the filtered data for 2018
cook_census_geo_2018_filtered_clipped.plot(column='B25077_001E', ax=ax[1], legend=True, scheme='Quantiles', cmap='coolwarm')
ax[1].set_title('Median Value of Owner-Occupied Housing Units (2018)', fontsize=16)
# Show the plot
plt.show()

from libpysal.weights import KNN
from esda.moran import Moran_Local

w_2018 = KNN.from_dataframe(cook_census_geo_2018_filtered_clipped, k=8)
w_2018.transform = "R"
lisa_2018 = Moran_Local(cook_census_geo_2018_filtered_clipped["B25077_001E"], w_2018)

w_2021 = KNN.from_dataframe(cook_census_geo_2021_filtered_clipped, k=8)
w_2021.transform = "R"
lisa_2021 = Moran_Local(cook_census_geo_2021_filtered_clipped["B25077_001E"], w_2021)

fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(24, 12))

cook_census_geo_2018_filtered_clipped.assign(Is=lisa_2018.Is).plot(
    column='Is',
    cmap='plasma',
    scheme='quantiles',
    k=5,
    edgecolor='white',
    linewidth=0.1,
    alpha=0.75,
    legend=True,
    ax=axs[0]
)
axs[0].set_title('Local Moran\'s I Values for Median Value of Owner-Occupied Housing Units - 2018', fontsize=16)
axs[0].set_axis_off()

cook_census_geo_2021_filtered_clipped.assign(Is=lisa_2021.Is).plot(
    column='Is',
    cmap='plasma',
    scheme='quantiles',
    k=5,
    edgecolor='white',
    linewidth=0.1,
    alpha=0.75,
    legend=True,
    ax=axs[1]
)
axs[1].set_title('Local Moran\'s I Values for Median Value of Owner-Occupied Housing Units - 2021', fontsize=16)
axs[1].set_axis_off()

# 显示整个图形
plt.show()

# Load 2023 CTA Bus Routes
cta_bus_routes = gpd.read_file(r"data\2023_CTA_Bus_Routes.geojson")
# Load 2023 CTA Bus Stops
cta_bus_stops = gpd.read_file(r"data\2023_CTA_Bus_Stops.geojson")
# Load 2023 CTA Rail Lines
cta_rail_lines = gpd.read_file(r"data\2023_CTA_Rail_Lines.geojson")
# Load 2023 CTA Rail Stations
cta_rail_stations = gpd.read_file(r"data\2023_CTA_Rail_Stations.geojson")
# Load 2023 Metra Rail Stations
metra_rail_stations = gpd.read_file(r"data\2023_Metra_Rail_Stations.geojson")
# Load 2023 Metra Rail Lines
metra_rail_lines = gpd.read_file(r"data\2023_Metra_Rail_Lines.geojson")
# Load 2023 Pace Bus Stops
pace_bus_stops = gpd.read_file(r"data\2023_Pace_Bus_Stops.geojson")
# Load 2023 Pace Bus Routes
pace_bus_routes = gpd.read_file(r"data\2023_Pace_Bus_Routes.geojson")
# Load Boundaries - Neighborhoods
neighborhood_boundaries = gpd.read_file(r"data\Boundaries - Neighborhoods.geojson")

cta_bus_routes.head(1)

cta_bus_stops.head(1)

# Create a figure and axis object
fig, ax = plt.subplots(figsize=(10, 10))
# Visualize boundary
cook_census_geo_2021_filtered_clipped.plot(ax=ax, color='#f5f5f5', edgecolor='grey', linewidth=0.5, label='Boundary')
# Visualize routes data (placed behind stops)
cta_bus_routes.plot(ax=ax, color='skyblue', linewidth=3,alpha=0.5,label='Routes')
# Visualize stops data (placed in front of routes)
cta_bus_stops.plot(ax=ax, color='blue', markersize=2, alpha=0.5,label='Stops')
# Add legend
ax.legend()
# Show the plot
plt.show()

# Create a figure and axis object
fig, ax = plt.subplots(figsize=(10, 10))
# Visualize boundary
cook_census_geo_2021_filtered_clipped.plot(ax=ax, color='#f5f5f5', edgecolor='grey', linewidth=0.5, label='Boundary')
# Visualize routes data (color by 'ROUTE' field)
cta_bus_routes.plot(ax=ax, column='ROUTE', linewidth=1, legend=True, zorder=1, label='Routes')
# Add legend
ax.legend()
# Show the plot
plt.show()

cta_rail_lines.head(1)

cta_rail_stations.head(1)

# Create a figure and axis object
fig, ax = plt.subplots(figsize=(10, 10))
# Visualize boundary
cook_census_geo_2021_filtered_clipped.plot(ax=ax, color='#f5f5f5', edgecolor='grey', linewidth=0.5, label='Boundary')
# Visualize routes data (placed behind stops)
cta_rail_lines.plot(ax=ax, color='skyblue', linewidth=3,alpha=0.5,label='Routes')
# Visualize stops data (placed in front of routes)
cta_rail_stations.plot(ax=ax, color='blue', markersize=4, alpha=0.5,label='Stops')
# Add legend
ax.legend()
# Show the plot
plt.show()

# Create a figure and axis object
fig, ax = plt.subplots(figsize=(10, 10))
# Visualize boundary
cook_census_geo_2021_filtered_clipped.plot(ax=ax, color='#f5f5f5', edgecolor='grey', linewidth=0.5, label='Boundary')
# Visualize routes data (color by 'ROUTE' field)
cta_rail_lines.plot(ax=ax, column='LINES', linewidth=3, legend=True, zorder=1, label='lines')
# Add legend
ax.legend()
# Show the plot
plt.show()

pace_bus_routes.head(1)

# Create a figure and axis object
fig, ax = plt.subplots(figsize=(10, 10))
# Visualize boundary
cook_census_geo_2021_filtered_clipped.plot(ax=ax, color='#f5f5f5', edgecolor='grey', linewidth=0.5, label='Boundary')
# Visualize routes data (placed behind stops)
pace_bus_routes.plot(ax=ax, color='skyblue', linewidth=3,alpha=0.5,label='Routes')
# Visualize stops data (placed in front of routes)
pace_bus_stops.plot(ax=ax, color='blue', markersize=2, alpha=0.5,label='Stops')
# Add legend
ax.legend()
# Show the plot
plt.show()

# Create a figure and axis object
fig, ax = plt.subplots(figsize=(10, 10))
# Visualize boundary
cook_census_geo_2021_filtered_clipped.plot(ax=ax, color='#f5f5f5', edgecolor='grey', linewidth=0.5, label='Boundary')
# Visualize routes data (color by 'ROUTE' field)
pace_bus_routes.plot(ax=ax, column='route_shor', linewidth=1, legend=True, zorder=1, label='Routes')
# Add legend
ax.legend()
# Show the plot
plt.show()

metra_rail_lines.head(1)

# Create a figure and axis object
fig, ax = plt.subplots(figsize=(10, 10))
# Visualize boundary
cook_census_geo_2021_filtered_clipped.plot(ax=ax, color='#f5f5f5', edgecolor='grey', linewidth=0.5, label='Boundary')
# Visualize routes data (placed behind stops)
metra_rail_lines.plot(ax=ax, color='skyblue', linewidth=3,alpha=0.5,label='Routes')
# Visualize stops data (placed in front of routes)
metra_rail_stations.plot(ax=ax, color='blue', markersize=4, alpha=0.5,label='Stops')
# Add legend
ax.legend()
# Show the plot
plt.show()

# Create a figure and axis object
fig, ax = plt.subplots(figsize=(10, 10))
# Visualize boundary
cook_census_geo_2021_filtered_clipped.plot(ax=ax, color='#f5f5f5', edgecolor='grey', linewidth=0.5, label='Boundary')
# Visualize routes data (color by 'ROUTE' field)
metra_rail_lines.plot(ax=ax, column='FID', linewidth=3, legend=True, zorder=1, label='lines')
# Add legend
ax.legend()
# Show the plot
plt.show()

file_path = r"data\CTA_Average_Rail_Station_Ridership_1999_2023.csv"
cta_ridership = pd.read_csv(file_path)
cta_ridership_2018 = cta_ridership[(cta_ridership['YEAR'] == 2018) & (cta_ridership['DAY_TYPE'] == 'Weekday')]
cta_ridership_2018 = cta_ridership_2018.groupby('RIDERSHIP_ID').agg({'DAILY_AVG_RIDES': 'mean'})
cta_ridership_2018.reset_index(inplace=True)
cta_ridership_2018.rename(columns={'RIDERSHIP_ID': 'STATION_ID'}, inplace=True)
cta_ridership_2018.head()

file_path = r"data\CTA_Average_Rail_Station_Ridership_1999_2023.csv"
cta_ridership = pd.read_csv(file_path)
cta_ridership_2021 = cta_ridership[(cta_ridership['YEAR'] == 2021) & (cta_ridership['DAY_TYPE'] == 'Weekday')]
cta_ridership_2021 = cta_ridership_2021.groupby('RIDERSHIP_ID').agg({'DAILY_AVG_RIDES': 'mean'})
cta_ridership_2021.reset_index(inplace=True)
cta_ridership_2021.rename(columns={'RIDERSHIP_ID': 'STATION_ID'}, inplace=True)
cta_ridership_2021.head()

file_path = r"data\Metra_Ridership_by_Station_Boarding_Alighting Survey_1979_2018_0.csv"
metra_ridership = pd.read_csv(file_path)
metra_ridership_2018 = metra_ridership[(metra_ridership['YEAR'] == 2018)]
metra_ridership_2018 = metra_ridership_2018.groupby('STATION').agg({'BOARDS': 'mean','ALIGHTS': 'mean'})
metra_ridership_2018.head()

cta_rail_stations.head(1)

cta_rail_stations.rename(columns={'LONGNAME': 'NAME'}, inplace=True)

cta_rail_stations.head(1)

cta_rail_ridership = pd.merge(cta_rail_stations, cta_ridership_2018[['STATION_ID', 'DAILY_AVG_RIDES']], left_on='STATION_ID', right_on='STATION_ID', how='left')
cta_rail_ridership.rename(columns={'DAILY_AVG_RIDES': 'RIDE2018'}, inplace=True)
cta_rail_ridership = pd.merge(cta_rail_ridership, cta_ridership_2021[['STATION_ID', 'DAILY_AVG_RIDES']], left_on='STATION_ID', right_on='STATION_ID', how='left')
cta_rail_ridership.rename(columns={'DAILY_AVG_RIDES': 'RIDE2021'}, inplace=True)
cta_rail_ridership.head()

cta_rail_ridership['railscale'] = cta_rail_ridership['RIDE2021'] / cta_rail_ridership['RIDE2018']
cta_rail_ridership.head()

cta_bus_routes.head(1)

#cta_bus_routes['ROUTE'].unique()

#cta_bus_routes['ROUTE0'].unique()

file_path = r"data\CTA_Average_Bus_Ridership_1999_2023.csv"
bus_ridership = pd.read_csv(file_path)
bus_ridership['ROUTE'].unique()

array(['0', '1', '100', '1001', '103', '106', '108', '11', '111', '111A',
       '112', '115', '119', '12', '120', '121', '124', '125', '126',
       '128', '134', '135', '136', '143', '146', '147', '148', '15',
       '151', '152', '155', '156', '157', '165', '169', '171', '172',
       '18', '19', '192', '2', '20', '201', '206', '21', '22', '24', '26',
       '28', '29', '3', '30', '31', '34', '35', '36', '37', '39', '4',
       '43', '44', '47', '48', '49', '49B', '5', '50', '51', '52', '52A',
       '53', '53A', '54', '54A', '54B', '55', '55A', '55N', '56', '57',
       '59', '6', '60', '62', '62H', '63', '63W', '65', '66', '67', '68',
       '7', '70', '71', '72', '73', '74', '75', '76', '77', '78', '79',
       '8', '80', '81', '81W', '82', '84', '85', '85A', '86', '87', '88',
       '8A', '9', '90', '91', '92', '93', '94', '95', '96', '97', 'J14',
       'X4', 'X49', 'X9', '10', '130', 'X98', '642', '132', '205', '95E',
       '95W', '170', '154', '33', 'R22', 'R39', 'R55', 'R63', 'R69',
       'R79', 'R87', 'R95', '122', '123', '129', '144', '145', '17',
       '49A', '56A', '64', '69', '90N', 'N201', 'X28', '14', '168',
       '53AL', 'X20', 'X3', 'X54', 'X55', 'X80', '174', '200', '38',
       '173', '127', '1002', '25', '40', 'X99', 'X21', '27', '202', '203',
       '204', 'N5', '1003'], dtype=object)

bus_ridership.head()

# Filter the data for the years 2018 and 2021
ridership_2018 = bus_ridership[bus_ridership['YEAR'] == 2018]
ridership_2021 = bus_ridership[bus_ridership['YEAR'] == 2021]

# Calculate the average monthly ridership by day type for each route
avg_ridership_2018 = ridership_2018.groupby(['ROUTE', 'DAY_TYPE'])['AVG_RIDES'].mean().unstack()
avg_ridership_2021 = ridership_2021.groupby(['ROUTE', 'DAY_TYPE'])['AVG_RIDES'].mean().unstack()

# Rename the columns as per your requirement
avg_ridership_2018.columns = ['saturday2018', 'sunday2018', 'weekday2018']
avg_ridership_2021.columns = ['saturday2021', 'sunday2021', 'weekday2021']

avg_ridership_2018.head()

avg_ridership_2021.head()

cta_bus_routes = pd.merge(cta_bus_routes, avg_ridership_2018, left_on='ROUTE', right_index=True, how='left')
cta_bus_routes = pd.merge(cta_bus_routes, avg_ridership_2021, left_on='ROUTE', right_index=True, how='left')

cta_bus_routes.head()

# Replace NaN values with 0
cta_bus_routes[['saturday2018', 'sunday2018', 'weekday2018', 'saturday2021', 'sunday2021', 'weekday2021']] = cta_bus_routes[['saturday2018', 'sunday2018', 'weekday2018', 'saturday2021', 'sunday2021', 'weekday2021']].fillna(0)
cta_bus_routes['weekend2018'] = cta_bus_routes['saturday2018'] + cta_bus_routes['sunday2018'] 
cta_bus_routes['weekend2021'] = cta_bus_routes['saturday2021'] + cta_bus_routes['sunday2021'] 
# Calculate total ridership for each year
cta_bus_routes['total2018'] = cta_bus_routes['saturday2018'] + cta_bus_routes['sunday2018'] + cta_bus_routes['weekday2018']
cta_bus_routes['total2021'] = cta_bus_routes['saturday2021'] + cta_bus_routes['sunday2021'] + cta_bus_routes['weekday2021']

# Calculate the scale
cta_bus_routes['busscale'] = cta_bus_routes['total2021'] / cta_bus_routes['total2018']
cta_bus_routes.head()

# Replace 0 with 1 in the specified columns of acs_df_2018 and acs_df_2021
acs_df_2018['B02001_001E'] = acs_df_2018['B02001_001E'].replace(0, 1)
acs_df_2021['B02001_001E'] = acs_df_2021['B02001_001E'].replace(0, 1)
acs_df_2018['B25077_001E'] = acs_df_2018['B25077_001E'].replace(0, 1)
acs_df_2021['B25077_001E'] = acs_df_2021['B25077_001E'].replace(0, 1)
acs_df_2018['B25085_001E'] = acs_df_2018['B25085_001E'].replace(0, 1)
acs_df_2021['B25085_001E'] = acs_df_2021['B25085_001E'].replace(0, 1)

# Calculate the logarithm of total population           
acs_df_2018['logPop'] = np.log(np.nan_to_num(acs_df_2018['B02001_001E']))
acs_df_2021['logPop'] = np.log(np.nan_to_num(acs_df_2021['B02001_001E']))

# Calculate log of housing values
acs_df_2018['logMidHomeValue'] = np.log(np.nan_to_num(acs_df_2018['B25077_001E']))
acs_df_2021['logMidHomeValue'] = np.log(np.nan_to_num(acs_df_2021['B25077_001E']))
acs_df_2018['logSaleHouseValue'] = np.log(np.nan_to_num(acs_df_2018['B25085_001E']))
acs_df_2021['logSaleHouseValue'] = np.log(np.nan_to_num(acs_df_2021['B25085_001E']))

acs_df_2018.head(1)

acs_df_2021.head(1)

print(acs_df_2018.columns)

Index(['NAME', 'B25077_001E', 'B25085_001E', 'B25001_001E', 'B25002_003E',
       'B25003_002E', 'B25003_003E', 'B06011_001E', 'B05007_001E',
       'B02001_001E', 'B02001_002E', 'state', 'county', 'tract', 'GEOID',
       'nonwhite_perc', 'foreign_perc', 'vacant_unit_perc',
       'owner_occupied_perc', 'renter_occupied_perc', 'logPop',
       'logMidHomeValue', 'logSaleHouseValue'],
      dtype='object')

acs_df_2018 = cook_tract.merge(acs_df_2018, left_on = 'GEOID', right_on = 'GEOID')
acs_df_2018.head()

print(acs_df_2018.columns)

Index(['STATEFP', 'COUNTYFP', 'TRACTCE', 'GEOID', 'NAME_x', 'NAMELSAD',
       'MTFCC', 'FUNCSTAT', 'ALAND', 'AWATER', 'INTPTLAT', 'INTPTLON',
       'geometry', 'NAME_y', 'B25077_001E', 'B25085_001E', 'B25001_001E',
       'B25002_003E', 'B25003_002E', 'B25003_003E', 'B06011_001E',
       'B05007_001E', 'B02001_001E', 'B02001_002E', 'state', 'county', 'tract',
       'nonwhite_perc', 'foreign_perc', 'vacant_unit_perc',
       'owner_occupied_perc', 'renter_occupied_perc', 'logPop',
       'logMidHomeValue', 'logSaleHouseValue'],
      dtype='object')

acs_df_2021 = cook_tract.merge(acs_df_2021, left_on = 'GEOID', right_on = 'GEOID')
acs_df_2021.head()

acs_df_2018.crs

<Geographic 2D CRS: EPSG:4269>
Name: NAD83
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: North America - onshore and offshore: Canada - Alberta; British Columbia; Manitoba; New Brunswick; Newfoundland and Labrador; Northwest Territories; Nova Scotia; Nunavut; Ontario; Prince Edward Island; Quebec; Saskatchewan; Yukon. Puerto Rico. United States (USA) - Alabama; Alaska; Arizona; Arkansas; California; Colorado; Connecticut; Delaware; Florida; Georgia; Hawaii; Idaho; Illinois; Indiana; Iowa; Kansas; Kentucky; Louisiana; Maine; Maryland; Massachusetts; Michigan; Minnesota; Mississippi; Missouri; Montana; Nebraska; Nevada; New Hampshire; New Jersey; New Mexico; New York; North Carolina; North Dakota; Ohio; Oklahoma; Oregon; Pennsylvania; Rhode Island; South Carolina; South Dakota; Tennessee; Texas; Utah; Vermont; Virginia; Washington; West Virginia; Wisconsin; Wyoming. US Virgin Islands. British Virgin Islands.
- bounds: (167.65, 14.92, -40.73, 86.45)
Datum: North American Datum 1983
- Ellipsoid: GRS 1980
- Prime Meridian: Greenwich

places = {
    'UChicago': [-87.60093, 41.78958],
    'UIC': [-87.65100, 41.87366],
    'Northwestern': [-87.67708, 42.05437],
    'NEIU': [-87.72031, 41.97870],
    'Illinois Tech': [-87.62721, 41.83139],
    'DePaul_Loop': [-87.62723, 41.87775],
    'DePaul_Lincoln_Park': [-87.65878, 41.92500],
    'Loyola_Lake_Shore': [-87.65779, 41.99858],
    'Loyola_Health_Sciences': [-87.83338, 41.86515],
    'Loyola_Water_Tower': [-87.62517, 41.89733],
    'Chicago_State': [-87.61350, 41.73000]
}

places_df = gpd.GeoDataFrame(
    {'place': list(places.keys()),
     'geometry': [Point(lon, lat) for lon, lat in places.values()]},
    crs="EPSG:4326"  # 使用WGS 84坐标系
)

acs_df_2018 = acs_df_2018.to_crs(epsg=3857)
places_df = places_df.to_crs(epsg=3857)

for place in places:
    acs_df_2018[f'distance_to_{place}'] = acs_df_2018.geometry.apply(
        lambda x: x.distance(places_df.loc[places_df['place'] == place, 'geometry'].squeeze())
    )

acs_df_2021 = acs_df_2021.to_crs(epsg=3857)
places_df = places_df.to_crs(epsg=3857)

for place in places:
    acs_df_2021[f'distance_to_{place}'] = acs_df_2021.geometry.apply(
        lambda x: x.distance(places_df.loc[places_df['place'] == place, 'geometry'].squeeze())
    )

acs_df_2018[['GEOID', 'INTPTLAT', 'INTPTLON', 'logPop',
       'logMidHomeValue', 'logSaleHouseValue', 'distance_to_UChicago',
       'distance_to_UIC', 'distance_to_Northwestern', 'distance_to_NEIU',
       'distance_to_Illinois Tech', 'distance_to_DePaul_Loop',
       'distance_to_DePaul_Lincoln_Park', 'distance_to_Loyola_Lake_Shore',
       'distance_to_Loyola_Health_Sciences', 'distance_to_Loyola_Water_Tower',
       'distance_to_Chicago_State']].head()

distances = np.array(acs_df_2018[['distance_to_UChicago', 'distance_to_UIC', 'distance_to_Northwestern', 'distance_to_NEIU', 'distance_to_Illinois Tech', 'distance_to_DePaul_Loop', 'distance_to_DePaul_Lincoln_Park', 'distance_to_Loyola_Lake_Shore', 'distance_to_Loyola_Health_Sciences', 'distance_to_Loyola_Water_Tower', 'distance_to_Chicago_State']])
acs_df_2018['dis_to_Uni'] = np.min(distances)
def get_min_distance(row):
    distances = [
        row['distance_to_UChicago'],
        row['distance_to_UIC'],
        row['distance_to_Northwestern'],
        row['distance_to_NEIU'],
        row['distance_to_Illinois Tech'],
        row['distance_to_DePaul_Loop'],
        row['distance_to_DePaul_Lincoln_Park'],
        row['distance_to_Loyola_Lake_Shore'],
        row['distance_to_Loyola_Health_Sciences'],
        row['distance_to_Loyola_Water_Tower'],
        row['distance_to_Chicago_State']
    ]
    return np.min(distances)

acs_df_2018['dis_to_Uni'] = acs_df_2018.apply(get_min_distance, axis=1)
# Create a list of the university distance columns
uni_distance_columns = ['distance_to_UChicago', 'distance_to_UIC', 'distance_to_Northwestern', 'distance_to_NEIU',
                        'distance_to_Illinois Tech', 'distance_to_DePaul_Loop', 'distance_to_DePaul_Lincoln_Park',
                        'distance_to_Loyola_Lake_Shore', 'distance_to_Loyola_Health_Sciences',
                        'distance_to_Loyola_Water_Tower', 'distance_to_Chicago_State']

# Create the 'closestUni' column by finding the column name with the minimum distance
acs_df_2018['closestUni'] = acs_df_2018[uni_distance_columns].idxmin(axis=1)
acs_df_2018['closestUni'] = acs_df_2018['closestUni'].str.replace('distance_to_', '')

distances = np.array(acs_df_2021[['distance_to_UChicago', 'distance_to_UIC', 'distance_to_Northwestern', 'distance_to_NEIU', 'distance_to_Illinois Tech', 'distance_to_DePaul_Loop', 'distance_to_DePaul_Lincoln_Park', 'distance_to_Loyola_Lake_Shore', 'distance_to_Loyola_Health_Sciences', 'distance_to_Loyola_Water_Tower', 'distance_to_Chicago_State']])
acs_df_2021['dis_to_Uni'] = np.min(distances)
def get_min_distance(row):
    distances = [
        row['distance_to_UChicago'],
        row['distance_to_UIC'],
        row['distance_to_Northwestern'],
        row['distance_to_NEIU'],
        row['distance_to_Illinois Tech'],
        row['distance_to_DePaul_Loop'],
        row['distance_to_DePaul_Lincoln_Park'],
        row['distance_to_Loyola_Lake_Shore'],
        row['distance_to_Loyola_Health_Sciences'],
        row['distance_to_Loyola_Water_Tower'],
        row['distance_to_Chicago_State']
    ]
    return np.min(distances)

acs_df_2021['dis_to_Uni'] = acs_df_2021.apply(get_min_distance, axis=1)
# Create a list of the university distance columns
uni_distance_columns = ['distance_to_UChicago', 'distance_to_UIC', 'distance_to_Northwestern', 'distance_to_NEIU',
                        'distance_to_Illinois Tech', 'distance_to_DePaul_Loop', 'distance_to_DePaul_Lincoln_Park',
                        'distance_to_Loyola_Lake_Shore', 'distance_to_Loyola_Health_Sciences',
                        'distance_to_Loyola_Water_Tower', 'distance_to_Chicago_State']

# Create the 'closestUni' column by finding the column name with the minimum distance
acs_df_2021['closestUni'] = acs_df_2021[uni_distance_columns].idxmin(axis=1)
acs_df_2021['closestUni'] = acs_df_2021['closestUni'].str.replace('distance_to_', '')

acs_df_2021.head(1)

private_campuses = ['UChicago', 'Northwestern', 'Illinois Tech', 'Loyola_Lake_Shore', 'Loyola_Health_Sciences', 'Loyola_Water_Tower']
acs_df_2018['is_private'] = acs_df_2018['closestUni'].apply(lambda x: 1 if x in private_campuses else 0)

acs_df_2021['is_private'] = acs_df_2021['closestUni'].apply(lambda x: 1 if x in private_campuses else 0)
acs_df_2021.head(1)

acs_df_2018['dis_to_Uni'].describe()

count     1288.000000
mean     10174.824797
std      10380.696481
min          0.000000
25%       3153.867030
50%       6482.344858
75%      12651.090871
max      50935.555851
Name: dis_to_Uni, dtype: float64

# Reproject cta_bus_stops and cta_rail_stations to a projected CRS
cta_bus_stops = cta_bus_stops.to_crs(epsg=3857)
cta_rail_stations = cta_rail_stations.to_crs(epsg=3857)
from shapely.ops import nearest_points

# Function to calculate the distance to the nearest bus stop
def calculate_distance_to_nearest_stop(row, stops_gdf):
    # Find the nearest bus stop
    nearest_stop = nearest_points(row['geometry'], stops_gdf.unary_union)[1]
    # Calculate the distance
    distance = row['geometry'].distance(nearest_stop)
    return distance

# Apply the function to each row in the acs_df GeoDataFrame
acs_df_2018['dis_to_bus'] = acs_df_2018.apply(calculate_distance_to_nearest_stop, axis=1, stops_gdf=cta_bus_stops)
acs_df_2018['dis_to_rail'] = acs_df_2018.apply(calculate_distance_to_nearest_stop, axis=1, stops_gdf=cta_rail_stations)
acs_df_2021['dis_to_bus'] = acs_df_2021.apply(calculate_distance_to_nearest_stop, axis=1, stops_gdf=cta_bus_stops)
acs_df_2021['dis_to_rail'] = acs_df_2021.apply(calculate_distance_to_nearest_stop, axis=1, stops_gdf=cta_rail_stations)

# Function to calculate the distance to the nearest rail station and get its ID and ridership
def calculate_distance_and_ridership(row, stations_gdf, ridership_gdf):
    # Find the nearest rail station
    nearest_station = nearest_points(row['geometry'], stations_gdf.unary_union)[1]
    # Get the station ID
    station_id = stations_gdf[stations_gdf['geometry'] == nearest_station]['STATION_ID'].values[0]
    # Get the ridership for the corresponding station ID
    ridership_2018 = ridership_gdf[ridership_gdf['STATION_ID'] == station_id]['RIDE2018'].values[0]
    ridership_2021 = ridership_gdf[ridership_gdf['STATION_ID'] == station_id]['RIDE2021'].values[0]
    # Calculate the distance
    distance = row['geometry'].distance(nearest_station)
    return pd.Series([distance, station_id, ridership_2018, ridership_2021], index=['dis_to_rail', 'station_id', 'ridership_2018', 'ridership_2021'])

# Apply the function to each row in the acs_df GeoDataFrame
acs_df_2018[['dis_to_rail', 'station_id', 'ridership_2018', 'ridership_2021']] = acs_df_2018.apply(calculate_distance_and_ridership, axis=1, stations_gdf=cta_rail_stations, ridership_gdf=cta_rail_ridership)
acs_df_2021[['dis_to_rail', 'station_id', 'ridership_2018', 'ridership_2021']] = acs_df_2021.apply(calculate_distance_and_ridership, axis=1, stations_gdf=cta_rail_stations, ridership_gdf=cta_rail_ridership)

# Load the Central Business District GeoJSON file
cbd_geojson = gpd.read_file("data\Central_Business_District_20240517.geojson")

# Ensure both GeoDataFrames use the same CRS
cbd_geojson = cbd_geojson.to_crs(acs_df_2018.crs)
# Calculate the centroid of the CBD if it's a polygon
cbd_centroid = cbd_geojson.geometry.centroid.iloc[0]
# Calculate the distance from each tract to the CBD centroid
acs_df_2018['dis_to_CBD'] = acs_df_2018.geometry.distance(cbd_centroid)
# Dummy 'is_CBD'
joined_df = gpd.sjoin(acs_df_2018, cbd_geojson, how="left", op="intersects")
joined_df['is_CBD'] = joined_df['index_right'].apply(lambda x: 1 if pd.notnull(x) else 0)
acs_df_2018 = acs_df_2018.merge(joined_df[['GEOID', 'is_CBD']], on='GEOID', how='left')

#2021
acs_df_2021['dis_to_CBD'] = acs_df_2021.geometry.distance(cbd_centroid)
# Dummy 'is_CBD'
joined_df = gpd.sjoin(acs_df_2021, cbd_geojson, how="left", op="intersects")
joined_df['is_CBD'] = joined_df['index_right'].apply(lambda x: 1 if pd.notnull(x) else 0)
acs_df_2021 = acs_df_2021.merge(joined_df[['GEOID', 'is_CBD']], on='GEOID', how='left')

acs_df_2021.head(1)

from esda.moran import Moran
import libpysal

# Assuming acs_df_2018 and acs_df_2021 are your GeoDataFrames for the years 2018 and 2021

# List of columns to include in the Moran's plot matrix
columns = [
    'logPop', 'nonwhite_perc', 'foreign_perc', 'vacant_unit_perc', 
    'owner_occupied_perc', 'dis_to_Uni', 'dis_to_bus', 'dis_to_rail'
]

# Function to create a spatial weights matrix
def create_weights(df):
    return libpysal.weights.Queen.from_dataframe(df)

# Function to calculate Moran's I for all columns
def calculate_all_morans_i(df, columns, w):
    return {column: Moran(df[column], w).I for column in columns}

# Create spatial weights matrices for both years
w_2018 = create_weights(acs_df_2018)
w_2021 = create_weights(acs_df_2021)

# Calculate Moran's I for all columns for both years
morans_i_2018 = calculate_all_morans_i(acs_df_2018, columns, w_2018)
morans_i_2021 = calculate_all_morans_i(acs_df_2021, columns, w_2021)

# Create the subplot grid with 4x4
fig, axs = plt.subplots(4, 4, figsize=(30, 20))  # Adjusted the size for better visibility

# Loop to plot each variable
for i, col in enumerate(columns):
    row = i // 4
    col_idx = i % 4
    sns.regplot(x=col, y='logMidHomeValue', data=acs_df_2018, ax=axs[row*2, col_idx], scatter_kws={'s': 0.5})
    axs[row*2, col_idx].set_title(f"2018 {col} (Moran's I: {morans_i_2018[col]:.4f})", fontsize=14)
    axs[row*2, col_idx].tick_params(axis='both', which='major', labelsize=12)

    sns.regplot(x=col, y='logMidHomeValue', data=acs_df_2021, ax=axs[row*2 + 1, col_idx], scatter_kws={'s': 0.5})
    axs[row*2 + 1, col_idx].set_title(f"2021 {col} (Moran's I: {morans_i_2021[col]:.4f})", fontsize=14)
    axs[row*2 + 1, col_idx].tick_params(axis='both', which='major', labelsize=12)

# Adjust layout and add a figure-wide title
plt.tight_layout()

plt.show()

variable_names = ['nonwhite_perc', 'foreign_perc', 'vacant_unit_perc', 'dis_to_rail']

Varlist18= variable_names + ["ridership_2018"]
acs_df_2018[acs_df_2018[Varlist18].isnull().any(axis=1)]

acs_df_2018 = acs_df_2018.dropna(subset=Varlist18+["logMidHomeValue"])

correlation_matrix = acs_df_2018[Varlist18].corr()
plt.figure(figsize=(10, 8))
sns.heatmap(correlation_matrix, annot=True, fmt=".2f")
plt.show()

from statsmodels.stats.outliers_influence import variance_inflation_factor
vif_data = pd.DataFrame()
vif_data["feature"] = acs_df_2018[Varlist18].columns
vif_data["VIF"] = [variance_inflation_factor(acs_df_2018[Varlist18].values, i) for i in range(len(acs_df_2018[Varlist18].columns))]
vif_data

OLS18 = spreg.OLS(
    acs_df_2018[["logMidHomeValue"]].values,
    acs_df_2018[Varlist18].values,
    name_y="logMidHomeValue",
    name_x=Varlist18,
)
print(OLS18.summary)

REGRESSION RESULTS
------------------

SUMMARY OF OUTPUT: ORDINARY LEAST SQUARES
-----------------------------------------
Data set            :     unknown
Weights matrix      :        None
Dependent Variable  :logMidHomeValue                Number of Observations:        1273
Mean dependent var  :     12.3303                Number of Variables   :           6
S.D. dependent var  :      0.5604                Degrees of Freedom    :        1267
R-squared           :      0.5475
Adjusted R-squared  :      0.5457
Sum squared residual:     180.722                F-statistic           :    306.6374
Sigma-square        :       0.143                Prob(F-statistic)     :  3.197e-215
S.E. of regression  :       0.378                Log likelihood        :    -563.753
Sigma-square ML     :       0.142                Akaike info criterion :    1139.505
S.E of regression ML:      0.3768                Schwarz criterion     :    1170.400

------------------------------------------------------------------------------------
            Variable     Coefficient       Std.Error     t-Statistic     Probability
------------------------------------------------------------------------------------
            CONSTANT        13.40960         0.03815       351.51406         0.00000
       nonwhite_perc        -1.16660         0.04251       -27.44573         0.00000
        foreign_perc        -0.78300         0.07915        -9.89314         0.00000
    vacant_unit_perc        -1.38949         0.17113        -8.11926         0.00000
         dis_to_rail        -0.00002         0.00000        -9.92575         0.00000
      ridership_2018        -0.00003         0.00000        -6.67313         0.00000
------------------------------------------------------------------------------------

REGRESSION DIAGNOSTICS
MULTICOLLINEARITY CONDITION NUMBER            8.662

TEST ON NORMALITY OF ERRORS
TEST                             DF        VALUE           PROB
Jarque-Bera                       2         225.992           0.0000

DIAGNOSTICS FOR HETEROSKEDASTICITY
RANDOM COEFFICIENTS
TEST                             DF        VALUE           PROB
Breusch-Pagan test                5         134.065           0.0000
Koenker-Bassett test              5          66.112           0.0000
================================ END OF REPORT =====================================

Varlist21= variable_names + ["ridership_2021"]
acs_df_2021 = acs_df_2021.dropna(subset=Varlist21+["logMidHomeValue"])
acs_df_2021[acs_df_2021[Varlist21].isnull().any(axis=1)]

correlation_matrix = acs_df_2021[Varlist21].corr()
plt.figure(figsize=(10, 8))
sns.heatmap(correlation_matrix, annot=True, fmt=".2f")
plt.show()

vif_data["feature"] = acs_df_2021[Varlist21].columns
vif_data["VIF"] = [variance_inflation_factor(acs_df_2021[Varlist21].values, i) for i in range(len(acs_df_2021[Varlist21].columns))]
vif_data

OLS21 = spreg.OLS(
    acs_df_2021[["logMidHomeValue"]].values,
    acs_df_2021[Varlist21].values,
    name_y="logMidHomeValue",
    name_x=Varlist21,
)
print(OLS21.summary)

REGRESSION RESULTS
------------------

SUMMARY OF OUTPUT: ORDINARY LEAST SQUARES
-----------------------------------------
Data set            :     unknown
Weights matrix      :        None
Dependent Variable  :logMidHomeValue                Number of Observations:        1312
Mean dependent var  :     12.4445                Number of Variables   :           6
S.D. dependent var  :      0.5577                Degrees of Freedom    :        1306
R-squared           :      0.5546
Adjusted R-squared  :      0.5529
Sum squared residual:     181.573                F-statistic           :    325.2886
Sigma-square        :       0.139                Prob(F-statistic)     :  2.122e-226
S.E. of regression  :       0.373                Log likelihood        :    -564.309
Sigma-square ML     :       0.138                Akaike info criterion :    1140.619
S.E of regression ML:      0.3720                Schwarz criterion     :    1171.695

------------------------------------------------------------------------------------
            Variable     Coefficient       Std.Error     t-Statistic     Probability
------------------------------------------------------------------------------------
            CONSTANT        13.51249         0.03564       379.17654         0.00000
       nonwhite_perc        -1.25649         0.04149       -30.28275         0.00000
        foreign_perc        -0.47560         0.07738        -6.14611         0.00000
    vacant_unit_perc        -1.19184         0.17808        -6.69273         0.00000
         dis_to_rail        -0.00001         0.00000        -8.75023         0.00000
      ridership_2021        -0.00006         0.00001        -6.96717         0.00000
------------------------------------------------------------------------------------

REGRESSION DIAGNOSTICS
MULTICOLLINEARITY CONDITION NUMBER            8.319

TEST ON NORMALITY OF ERRORS
TEST                             DF        VALUE           PROB
Jarque-Bera                       2         671.549           0.0000

DIAGNOSTICS FOR HETEROSKEDASTICITY
RANDOM COEFFICIENTS
TEST                             DF        VALUE           PROB
Breusch-Pagan test                5          68.764           0.0000
Koenker-Bassett test              5          25.558           0.0001
================================ END OF REPORT =====================================

SR18 = spreg.OLS_Regimes(
    acs_df_2018[["logMidHomeValue"]].values,
    acs_df_2018[Varlist18].values,
    acs_df_2018["is_CBD"].tolist(),
    constant_regi="many",
    regime_err_sep=False,
    name_y="logMidHomeValue",
    name_x=Varlist18,
)

res = pd.DataFrame(
    {
        # Pull out regression coefficients and
        # flatten as they are returned as Nx1 array
        "Coeff.": SR18.betas.flatten(),
        # Pull out and flatten standard errors
        "Std. Error": SR18.std_err.flatten(),
        # Pull out P-values from t-stat object
        "P-Value": [i[1] for i in SR18.t_stat],
    },
    index=SR18.name_x,
)
## Extract variables for the coastal regime
CBD = [i for i in res.index if "1_" in i]
CBD = res.loc[CBD, :].rename(lambda i: i.replace("1_", ""))
## Build multi-index column names
CBD.columns = pd.MultiIndex.from_product(
    [["CBD"], CBD.columns]
)
# Non-CBD
## Extract variables for the non-coastal regime
out = [i for i in res.index if "0_" in i]
out = res.loc[out, :].rename(lambda i: i.replace("0_", ""))
## Build multi-index column names
out.columns = pd.MultiIndex.from_product(
    [["Outside"], out.columns]
)
# Concat both models
pd.concat([CBD, out], axis=1)

SR21 = spreg.OLS_Regimes(
    acs_df_2021[["logMidHomeValue"]].values,
    acs_df_2021[Varlist21].values,
    acs_df_2021["is_CBD"].tolist(),
    constant_regi="many",
    regime_err_sep=False,
    name_y="logMidHomeValue",
    name_x=Varlist21,
)

res = pd.DataFrame(
    {
        # Pull out regression coefficients and
        # flatten as they are returned as Nx1 array
        "Coeff.": SR21.betas.flatten(),
        # Pull out and flatten standard errors
        "Std. Error": SR21.std_err.flatten(),
        # Pull out P-values from t-stat object
        "P-Value": [i[1] for i in SR21.t_stat],
    },
    index=SR21.name_x,
)
## Extract variables for the coastal regime
CBD = [i for i in res.index if "1_" in i]
CBD = res.loc[CBD, :].rename(lambda i: i.replace("1_", ""))
## Build multi-index column names
CBD.columns = pd.MultiIndex.from_product(
    [["CBD"], CBD.columns]
)
# Non-CBD
## Extract variables for the non-coastal regime
out = [i for i in res.index if "0_" in i]
out = res.loc[out, :].rename(lambda i: i.replace("0_", ""))
## Build multi-index column names
out.columns = pd.MultiIndex.from_product(
    [["Outside"], out.columns]
)
# Concat both models
pd.concat([CBD, out], axis=1)

knn18 = weights.KNN.from_dataframe(acs_df_2018, k=20)
SLM18 = spreg.GM_Lag(
    acs_df_2018[["logMidHomeValue"]].values,
    acs_df_2018[Varlist18].values,
    w=knn18,
    name_y="logMidHomeValue",
    name_x=Varlist18,
)
print(SLM18.summary)

REGRESSION RESULTS
------------------

SUMMARY OF OUTPUT: SPATIAL TWO STAGE LEAST SQUARES
--------------------------------------------------
Data set            :     unknown
Weights matrix      :     unknown
Dependent Variable  :logMidHomeValue                Number of Observations:        1273
Mean dependent var  :     12.3303                Number of Variables   :           7
S.D. dependent var  :      0.5604                Degrees of Freedom    :        1266
Pseudo R-squared    :      0.7732
Spatial Pseudo R-squared:  0.6118

------------------------------------------------------------------------------------
            Variable     Coefficient       Std.Error     z-Statistic     Probability
------------------------------------------------------------------------------------
            CONSTANT         4.31300         0.51284         8.41001         0.00000
       nonwhite_perc        -0.49124         0.04854       -10.12075         0.00000
        foreign_perc        -0.61779         0.05694       -10.84917         0.00000
    vacant_unit_perc        -0.81194         0.12575        -6.45677         0.00000
         dis_to_rail        -0.00000         0.00000        -3.70018         0.00022
      ridership_2018        -0.00001         0.00000        -4.50443         0.00001
   W_logMidHomeValue         0.03468         0.00195        17.76240         0.00000
------------------------------------------------------------------------------------
Instrumented: W_logMidHomeValue
Instruments: W_dis_to_rail, W_foreign_perc, W_nonwhite_perc,
             W_ridership_2018, W_vacant_unit_perc
================================ END OF REPORT =====================================

knn21 = weights.KNN.from_dataframe(acs_df_2021, k=20)
SLM21 = spreg.GM_Lag(
    acs_df_2021[["logMidHomeValue"]].values,
    acs_df_2021[Varlist21].values,
    w=knn21,
    name_y="logMidHomeValue",
    name_x=Varlist21,
)
print(SLM21.summary)

REGRESSION RESULTS
------------------

SUMMARY OF OUTPUT: SPATIAL TWO STAGE LEAST SQUARES
--------------------------------------------------
Data set            :     unknown
Weights matrix      :     unknown
Dependent Variable  :logMidHomeValue                Number of Observations:        1312
Mean dependent var  :     12.4445                Number of Variables   :           7
S.D. dependent var  :      0.5577                Degrees of Freedom    :        1305
Pseudo R-squared    :      0.7552
Spatial Pseudo R-squared:  0.6030

------------------------------------------------------------------------------------
            Variable     Coefficient       Std.Error     z-Statistic     Probability
------------------------------------------------------------------------------------
            CONSTANT         4.45694         0.56793         7.84765         0.00000
       nonwhite_perc        -0.52914         0.05501        -9.61847         0.00000
        foreign_perc        -0.48426         0.05749        -8.42402         0.00000
    vacant_unit_perc        -0.76462         0.13497        -5.66532         0.00000
         dis_to_rail        -0.00000         0.00000        -3.05028         0.00229
      ridership_2021        -0.00003         0.00001        -4.08048         0.00004
   W_logMidHomeValue         0.03421         0.00214        15.96210         0.00000
------------------------------------------------------------------------------------
Instrumented: W_logMidHomeValue
Instruments: W_dis_to_rail, W_foreign_perc, W_nonwhite_perc,
             W_ridership_2021, W_vacant_unit_perc
================================ END OF REPORT =====================================

knn18.transform = "R"
lag_variables_used = acs_df_2018.filter(['nonwhite_perc', 'foreign_perc', 'vacant_unit_perc', 'dis_to_rail', 'ridership_2018'])
wx = lag_variables_used.apply(lambda y: weights.spatial_lag.lag_spatial(knn18, y))
wx = wx.rename(columns=lambda c: "w_"+ c)
slx_exog18 = acs_df_2018[Varlist18].join(wx)

knn21.transform = "R"
lag_variables_used = acs_df_2021.filter(['nonwhite_perc', 'foreign_perc', 'vacant_unit_perc', 'dis_to_rail', 'ridership_2021'])
wx = lag_variables_used.apply(lambda y: weights.spatial_lag.lag_spatial(knn21, y))
wx = wx.rename(columns=lambda c: "w_"+ c)
slx_exog21 = acs_df_2021[Varlist21].join(wx)

SDM18 = spreg.GM_Lag(
    acs_df_2018[["logMidHomeValue"]].values,
    slx_exog18.values,
    w=knn18,
    name_y="logMidHomeValue",
    name_x=slx_exog18.columns.tolist(),
)
print(SDM18.summary)

REGRESSION RESULTS
------------------

SUMMARY OF OUTPUT: SPATIAL TWO STAGE LEAST SQUARES
--------------------------------------------------
Data set            :     unknown
Weights matrix      :     unknown
Dependent Variable  :logMidHomeValue                Number of Observations:        1273
Mean dependent var  :     12.3303                Number of Variables   :          12
S.D. dependent var  :      0.5604                Degrees of Freedom    :        1261
Pseudo R-squared    :      0.0202
Spatial Pseudo R-squared: omitted due to rho outside the boundary (-1, 1).

------------------------------------------------------------------------------------
            Variable     Coefficient       Std.Error     z-Statistic     Probability
------------------------------------------------------------------------------------
            CONSTANT       -34.42374    
       nonwhite_perc        -2.78033         5.08284        -0.54700         0.58438
        foreign_perc         2.70351         8.60099         0.31433         0.75327
    vacant_unit_perc        -2.36103        13.75619        -0.17163         0.86373
         dis_to_rail         0.00001         0.00063         0.01937         0.98454
      ridership_2018         0.00008         0.00052         0.14965         0.88104
     w_nonwhite_perc         6.50409    
      w_foreign_perc         1.11162        14.24899         0.07801         0.93782
  w_vacant_unit_perc         7.89387    
       w_dis_to_rail        -0.00002         0.00007        -0.23264         0.81604
    w_ridership_2018         0.00001         0.00020         0.03967         0.96835
   W_logMidHomeValue         1.45985    
------------------------------------------------------------------------------------
Instrumented: W_logMidHomeValue
Instruments: W_dis_to_rail, W_foreign_perc, W_nonwhite_perc,
             W_ridership_2018, W_vacant_unit_perc, W_w_dis_to_rail,
             W_w_foreign_perc, W_w_nonwhite_perc, W_w_ridership_2018,
             W_w_vacant_unit_perc
Warning: *** WARNING: Estimate for spatial lag coefficient is outside the boundary (-1, 1). ***
================================ END OF REPORT =====================================

	ALAND	AWATER	B25077_001E	B25085_001E	B25001_001E	B25002_003E	B25003_002E	B25003_003E	B06011_001E	B05007_001E	B02001_001E	B02001_002E	nonwhite_perc	foreign_perc	vacant_unit_perc	owner_occupied_perc	renter_occupied_perc
count	1.332000e+03	1.332000e+03	1.332000e+03	1332.000000	1332.000000	1332.000000	1332.000000	1332.000000	1.332000e+03	1332.000000	1332.000000	1332.000000	1328.000000	1328.000000	1328.000000	1328.000000	1328.000000
mean	1.837350e+06	1.341090e+06	-9.719136e+06	20.870120	1692.144144	157.115616	882.239489	652.789039	-3.463530e+06	827.337087	3953.001502	2057.081081	0.511296	0.192016	0.097898	0.515126	0.386976
std	3.147841e+06	4.704685e+07	8.114146e+07	28.637381	835.675693	138.700517	564.369966	533.515014	4.822262e+07	778.174197	1778.569024	1577.538675	0.294288	0.142000	0.071966	0.233497	0.201916
min	0.000000e+00	0.000000e+00	-6.666667e+08	0.000000	0.000000	0.000000	0.000000	0.000000	-6.666667e+08	0.000000	0.000000	0.000000	0.034380	0.000000	0.000000	0.000000	0.000000
25%	4.110900e+05	0.000000e+00	1.736750e+05	0.000000	1124.750000	70.000000	423.500000	295.000000	2.683000e+04	210.000000	2560.250000	656.000000	0.255031	0.071450	0.048706	0.326816	0.220682
50%	8.969725e+05	0.000000e+00	2.483000e+05	10.000000	1614.500000	122.000000	794.000000	526.000000	3.472350e+04	589.000000	3885.000000	1895.500000	0.453815	0.165461	0.081027	0.500264	0.386741
75%	2.147472e+06	0.000000e+00	3.627000e+05	33.000000	2114.000000	203.250000	1264.500000	833.000000	4.929225e+04	1256.000000	5145.750000	3170.250000	0.777141	0.300599	0.125492	0.711206	0.541907
max	6.184753e+07	1.717072e+09	1.474100e+06	219.000000	7694.000000	1345.000000	3108.000000	4750.000000	1.147500e+05	4543.000000	11295.000000	8515.000000	1.000000	0.669796	0.477829	1.000000	0.985468

	ALAND	AWATER	B25077_001E	B25085_001E	B25001_001E	B25002_003E	B25003_002E	B25003_003E	B06011_001E	B05007_001E	B02001_001E	B02001_002E	nonwhite_perc	foreign_perc	vacant_unit_perc	owner_occupied_perc	renter_occupied_perc
count	1.288000e+03	1.288000e+03	1.288000e+03	1288.000000	1288.000000	1288.000000	1288.000000	1288.000000	1.288000e+03	1288.000000	1288.000000	1288.000000	1284.000000	1284.000000	1284.000000	1284.000000	1284.000000
mean	1.742769e+06	1.382477e+06	-7.502018e+06	28.131988	1635.746118	169.021739	829.319099	637.405280	-2.036330e+06	821.369565	3914.440217	2212.215839	0.469626	0.191050	0.109513	0.496794	0.393693
std	2.648591e+06	4.784372e+07	7.158056e+07	33.609735	771.371540	138.654729	541.363786	478.766688	3.711046e+07	789.686418	1732.385134	1663.625347	0.314474	0.145032	0.078187	0.230792	0.194565
min	0.000000e+00	0.000000e+00	-6.666667e+08	0.000000	0.000000	0.000000	0.000000	0.000000	-6.666667e+08	0.000000	0.000000	0.000000	0.009518	0.000000	0.000000	0.000000	0.008403
25%	4.019952e+05	0.000000e+00	1.533000e+05	0.000000	1106.000000	76.000000	386.000000	306.000000	2.235900e+04	201.000000	2546.250000	685.750000	0.188690	0.067770	0.055406	0.314252	0.235574
50%	8.629060e+05	0.000000e+00	2.215500e+05	19.000000	1576.000000	136.000000	750.500000	525.500000	3.000000e+04	569.500000	3833.000000	2140.500000	0.398953	0.162402	0.090442	0.474632	0.400432
75%	2.038647e+06	0.000000e+00	3.323750e+05	43.000000	2058.000000	217.250000	1183.000000	829.500000	4.129100e+04	1266.000000	5138.250000	3397.500000	0.758785	0.305400	0.140072	0.687850	0.546442
max	4.514031e+07	1.717072e+09	1.528600e+06	311.000000	6514.000000	1140.000000	3382.000000	3686.000000	1.057860e+05	4893.000000	10317.000000	8467.000000	1.000000	0.691883	0.569182	0.971333	0.942693

	STATEFP	COUNTYFP	TRACTCE	GEOID	NAME_x	NAMELSAD	MTFCC	FUNCSTAT	ALAND	...	B02001_001E	B02001_002E	state	county	tract	nonwhite_perc	foreign_perc	vacant_unit_perc	owner_occupied_perc	renter_occupied_perc
0	17	031	221000	17031221000	2210	Census Tract 2210	G5020	S	220053	...	2774.0	2099.0	17	031	221000	0.243331	0.169430	0.089847	0.297394	0.612758
1	17	031	221100	17031221100	2211	Census Tract 2211	G5020	S	440042	...	5207.0	3900.0	17	031	221100	0.251008	0.288650	0.072526	0.286263	0.641210
2	17	031	242800	17031242800	2428	Census Tract 2428	G5020	S	439630	...	1542.0	1461.0	17	031	242800	0.052529	0.127756	0.079812	0.593114	0.327074
3	17	031	242900	17031242900	2429	Census Tract 2429	G5020	S	323819	...	1617.0	1229.0	17	031	242900	0.239951	0.136054	0.107987	0.388076	0.503937
4	17	031	243000	17031243000	2430	Census Tract 2430	G5020	S	324547	...	2274.0	1706.0	17	031	243000	0.249780	0.174582	0.040075	0.363467	0.596459

	BOARDS	ALIGHTS
STATION
103rd Street (Beverly Hills)	734.0	731.0
103rd Street (Rosemoor)	36.0	34.0
107th Street	27.0	23.0
107th Street (Beverly Hills)	395.0	368.0
111th Street	31.0	49.0

	FID	STATION_ID	NAME	LINES	ADDRESS	ADA	PKNRD	POINT_X	POINT_Y	geometry	RIDE2018	RIDE2021
0	1	230	Cumberland	Blue Line	5800 N. Cumberland Avenue	1	1	1.118914e+06	1.937256e+06	POINT (-87.83803 41.98430)	4498.583333	1177.333333
1	2	1350	Oak Park-Lake	Green Line (Lake)	100 S. Oak Park Avenue	0	0	1.131166e+06	1.901870e+06	POINT (-87.79379 41.88700)	1578.583333	426.166667
2	3	1260	Austin-Lake	Green Line (Lake)	351 N. Austin Blvd	0	0	1.136515e+06	1.902016e+06	POINT (-87.77414 41.88730)	1822.916667	701.583333
3	4	170	Ashland-Lake	Green (Lake), Pink	1601 W. Lake Street	1	0	1.165707e+06	1.901502e+06	POINT (-87.66696 41.88532)	2431.750000	992.583333
4	5	1160	Clinton-Lake	Green (Lake), Pink	540 W. Lake Street	1	0	1.172562e+06	1.901699e+06	POINT (-87.64178 41.88571)	4406.916667	1247.250000

	STATION_ID	DAILY_AVG_RIDES
0	10	1934.250000
1	20	3746.916667
2	30	1396.000000
3	40	7961.416667
4	50	3639.750000

	STATION_ID	DAILY_AVG_RIDES
0	10	551.250000
1	20	1404.833333
2	30	657.916667
3	40	1703.333333
4	50	1163.083333

	ROUTE	YEAR	MONTH	DAY_TYPE	AVG_RIDES
0	0	2023	10	Weekday	2626
1	0	2023	10	Saturday	4694
2	0	2023	10	Sunday - Holiday	2634
3	1	2023	10	Weekday	1126
4	100	2023	10	Weekday	470

	saturday2018	sunday2018	weekday2018
ROUTE
0	2710.833333	1586.50	3680.833333
1	NaN	NaN	1709.833333
10	455.600000	473.60	521.000000
100	NaN	NaN	569.666667
1001	139.083333	121.75	85.916667

	saturday2021	sunday2021	weekday2021
ROUTE
0	2079.666667	1540.250000	3063.500000
1	NaN	NaN	376.833333
10	408.000000	331.400000	362.250000
100	NaN	NaN	243.416667
1001	463.750000	495.583333	53.416667

	FID	ROUTE	ROUTE0	NAME	WKDAY	SAT	SUN	SHAPE_LEN	Shape__Length	geometry	saturday2018	sunday2018	weekday2018	saturday2021	sunday2021	weekday2021
0	1	12	012	ROOSEVELT	1	1	1	62586.111408	25595.796952	MULTILINESTRING ((-87.64792 41.86714, -87.6474...	7589.083333	5757.416667	12381.416667	3761.750000	2735.916667	5725.083333
1	2	121	121	UNION/STREETERVILLE EXPRESS	1	0	0	24090.823988	9866.558032	MULTILINESTRING ((-87.62451 41.88834, -87.6249...	NaN	NaN	1111.666667	NaN	NaN	130.250000
2	3	1	001	BRONZEVILLE/UNION STATION	1	0	0	34690.953676	14214.349988	MULTILINESTRING ((-87.62327 41.83104, -87.6232...	NaN	NaN	1709.833333	NaN	NaN	376.833333
3	4	108	108	HALSTED/95TH	1	0	0	71888.990825	29369.622395	MULTILINESTRING ((-87.59054 41.65578, -87.5906...	NaN	NaN	1138.500000	NaN	NaN	380.583333
4	5	11	011	LINCOLN	1	1	1	24694.573889	10138.959924	MULTILINESTRING ((-87.68884 41.96677, -87.6888...	799.583333	528.000000	1452.333333	491.416667	337.166667	767.000000

	GEOID	INTPTLAT	INTPTLON	logPop	logMidHomeValue	logSaleHouseValue	distance_to_UChicago	distance_to_UIC	distance_to_Northwestern	distance_to_NEIU	distance_to_Illinois Tech	distance_to_DePaul_Loop	distance_to_DePaul_Lincoln_Park	distance_to_Loyola_Lake_Shore	distance_to_Loyola_Health_Sciences	distance_to_Loyola_Water_Tower	distance_to_Chicago_State
0	17031221000	+41.9209604	-087.7152526	7.928046	12.570716	2.302585	22837.922146	9543.405442	19853.038273	8105.724679	16042.463991	11282.666455	6117.381881	12700.475695	15141.922884	10279.749562	30116.492337
1	17031221100	+41.9210205	-087.7103609	8.557759	12.746069	3.044522	22456.389068	9034.397115	19706.054739	8124.189957	15626.730254	10675.770660	5391.661358	12351.155543	15456.313556	9588.712436	29865.883552
2	17031242800	+41.8926320	-087.6903237	7.340836	13.010090	0.000000	17580.179771	4537.583407	23777.110734	12702.440882	10787.868719	6801.079518	5376.957816	15727.818981	15831.552924	6869.467102	25019.487371
3	17031242900	+41.8921955	-087.6842896	7.388328	12.937479	0.000000	17320.705974	4089.832277	23745.884872	12961.290725	10491.517010	6282.694874	5069.011206	15612.515897	16706.620747	6324.563528	24880.809413
4	17031243000	+41.8922753	-087.6793981	7.729296	12.892195	2.397895	17054.379038	3653.502848	23726.987406	13116.382128	10195.929179	5761.598949	4803.055148	15515.312166	17242.005834	5779.328112	24731.243256

	STATEFP	COUNTYFP	TRACTCE	GEOID	NAME_x	NAMELSAD	MTFCC	FUNCSTAT	ALAND	AWATER	...	dis_to_Uni	closestUni	is_private	dis_to_bus	dis_to_rail	station_id	ridership_2018	ridership_2021	dis_to_CBD	is_CBD
55	17	031	980100	17031980100	9801	Census Tract 9801	G5020	S	2981781	0	...	13402.991508	Loyola_Health_Sciences	1	0.000000	0.000000	930.0	8599.833333	3180.250000	18247.821379	0
323	17	031	980000	17031980000	9800	Census Tract 9800	G5020	S	19890200	92402	...	14399.445316	Loyola_Health_Sciences	1	3325.041188	0.000000	890.0	11400.916667	5028.833333	30080.554490	0
750	17	031	381700	17031381700	3817	Census Tract 3817	G5020	S	196802	0	...	3311.314365	Illinois Tech	1	0.000000	292.471936	1230.0	2879.500000	1110.916667	11248.677147	0
1161	17	031	990000	17031990000	9900	Census Tract 9900	G5020	S	0	1717072182	...	755.516274	Loyola_Lake_Shore	1	296.765334	1028.539570	1300.0	5129.833333	1953.583333	2233.874727	1

	feature	VIF
0	nonwhite_perc	4.199090
1	foreign_perc	1.908475
2	vacant_unit_perc	4.076822
3	dis_to_rail	2.522510
4	ridership_2018	4.002518

	feature	VIF
0	nonwhite_perc	4.442830
1	foreign_perc	2.061458
2	vacant_unit_perc	3.761674
3	dis_to_rail	2.854869
4	ridership_2021	4.462365

	CBD			Outside
	Coeff.	Std. Error	P-Value	Coeff.	Std. Error	P-Value
CONSTANT	12.972485	0.269899	2.508697e-287	13.455864	0.037573	0.000000e+00
nonwhite_perc	0.028624	0.843532	9.729353e-01	-1.100378	0.041823	5.689917e-122
foreign_perc	-1.023915	1.081883	3.441164e-01	-0.812831	0.076990	5.009319e-25
vacant_unit_perc	0.618738	0.838265	4.605805e-01	-1.851549	0.173497	1.627757e-25
dis_to_rail	0.000271	0.000317	3.934378e-01	-0.000013	0.000002	6.172307e-16
ridership_2018	0.000012	0.000023	6.145541e-01	-0.000035	0.000004	3.158383e-19

	CBD			Outside
	Coeff.	Std. Error	P-Value	Coeff.	Std. Error	P-Value
CONSTANT	13.266159	0.303604	2.258006e-257	13.528267	0.035486	0.000000e+00
nonwhite_perc	-0.906958	0.787932	2.499188e-01	-1.201198	0.041690	1.383110e-141
foreign_perc	0.110858	0.759801	8.840197e-01	-0.511979	0.076788	3.837076e-11
vacant_unit_perc	0.715346	0.988215	4.692716e-01	-1.525842	0.182153	1.393752e-16
dis_to_rail	0.000205	0.000305	5.016357e-01	-0.000012	0.000002	1.243038e-13
ridership_2021	-0.000029	0.000064	6.472906e-01	-0.000078	0.000009	2.536415e-16

Spatial Analysis of Demographic, Housing Prices and Public Transportation Accessibility in Chicago¶

Yizhou Wang & Hongtao Huang¶

Master of Regional Planning 25’¶

1 Introduction¶

2 Problem Statement¶

3 Literature Review¶

3.1 Housing Prices and Public Transportation¶

3.2 Public Transportation in Chicago¶

3.2.1 The Demographic of Public Transportation in Chicago¶

3.2.2 Influence of Covid to the Public Transit in Chicago¶

4 Methodology plan¶

4.1 Data:¶

4.1.1 Demographic and Economic data in Block Group: Census API¶

4.1.2 Public Transit Data in Chicago¶

4.2 Analysis methods¶

5 Site Analysis of Basic Data¶

5.1 Census Demographic Data¶

5.1 Data Processing¶

5.2 Visualization of Social-Economic Factors and Housing Price in Chicago¶

5.2.1 Social-Economic Pattern in Chicago¶

5.2.1 Housing Price Pattern in Chicago¶

5.2 Transit Routes and Stations in Chicago¶

5.2.1 CTA Bus Route and Stops¶

5.2.2 CTA Rail Lines and Stations¶

5.2.3 Pace Bus Route and Stops¶

5.2.4 Metra Rail and Lines¶

5.3 Data Merging of Ridership Data¶

5.3.1 CTA Ridership Data¶

5.3.2 Metra Ridership Data¶

Join ridership data to the subway station data¶

Join ridership data to the bus line data¶

Calculate CTA bus ridership by weekday, saturday, and sunday¶

Merge 2 years ridership data together¶

Calculate the weekend bus ridership¶

6 Models¶

6.1 Variables¶

6.2 Data Processing¶

6.3 Global Moran's I¶

6.3.1 Calculation for Global Moran's I¶

6.3.2 Interpretation of Global Moran's I Graphs¶

6.4 OLS Regression¶

6.4.1 OLS for 2018¶

6.4.2 OLS Model for 2021¶

6.4.3 Interpretation of OLS Models¶

6.5 Spatial Regime Model¶

6.5.1 Spatial Regime Model for 2018¶

6.5.2 Spatial Regime Model for 2021¶

6.5.3 Interpretation for Spatial Regime Model¶

6.6 Spatial Lag Model¶

6.6.1 Spatial Lag Model for 2018¶

6.6.2 Spatial Lag Model for 2021¶

6.6.3 Interpretation for Spatial Lag Model¶

6.6 Spatial Durbin Model¶

6.6.1 Spatial Durbin Model for 2018¶

6.6.2 Limitation for Spatial Durbin Model¶

7 Limitations¶

8 Conclusion¶

9 Reference¶